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5h ' We carried out global three-dimensional resistive magnetohydrodynamic sim- 

ulations of galactic gaseous disks to investigate how the galactic magnetic fields 
are amplified and maintained. We adopt a steady axisymmetric gravitational 
potential given by Miyamoto & Nagai (1975) and Miyamoto, Satoh & Ohashi 
(1980). As the initial condition, we assume a warm (T ~ 10 5 K) rotating gas 
torus centered at w = 10 kpc threaded by weak azimuthal magnetic fields. Nu- 
merical results indicate that in differentially rotating galactic gaseous disks, mag- 
netic fields are amplified due to magneto-rotational instability and magnetic tur- 
bulence develops. After the amplification of magnetic energy saturates, the disk 
stays in a quasi-steady state. The mean azimuthal magnetic field increases with 
time and shows reversals with period of lGyr (2Gyr for a full cycle). The am- 
plitude of B v near the equatorial plane is B v ~ 1.5 fi G at zu = 5 kpc. The 
magnetic fields show large fluctuations whose standard deviation is comparable 
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to the mean field. The mean azimuthal magnetic field in the disk corona has di- 
rection opposite to the mean magnetic field inside the disk. The mass accretion 
rate driven by the Maxwell stress is ~ 1O~ 3 M / yr at w = 2.5 kpc when the 
mass of the initial torus is ~ 5 x 1O 8 M . When we adopt an absorbing boundary 
condition at r = 0.8 kpc, the rotation curve obtained by numerical simulations 
almost coincides with the rotation curve of the stars and the dark matter. Thus 
even when magnetic fields are not negligible for gas dynamics of a spiral galaxy 
galactic gravitational potential can be derived from observation of rotation curves 
using gas component of the disk. 

Subject headings: galaxies: spiral — ISM: magnetic fields — MHD 



1. INTRODUCTION 

Spiral galaxies have magnetic fields. In our Galaxy the average magnetic field strength 
is a few fiG. (e.g., Sofue, Fujimoto & Wielebinski 1986; Beck et al. 1996). These magnetic 
fields have been explained by the dynamo action operating in galactic gas disks (e.g., Parker 
1971). In conventional theory of galactic dynamos, the following induction equation is solved; 

f) Pi 

— = V x (v x B) +r]V 2 B + V x (aB), (1) 

where 77 denotes the turbulent magnetic diffusivity, v is the mean velocity, and B is the mean 
magnetic field. The last term on the right hand side of equation (1) is the dynamo term which 
represents the re-generation of mean magnetic fields (a effect). When velocity fields v are 
given, equation (1) is linear with B. It can be solved with appropriate boundary conditions. 
This approach is referred to as the kinematic dynamo. The effects of nonlinearities such as 
the quenching of the ct-effect for strong magnetic fields and the loss of magnetic flux due to 
the magnetic buoyancy have also been studied in the framework of the kinematic dynamo 
theory (e.g., Schmitt & Schiissler 1989; Brandenburg et al. 1989, 1992). 

In 1990s, it became clear that nonlinear dynamical processes are essential for the evo- 
lution of magnetic fields in differentially rotating disks. Balbus & Hawley (1991) showed 
that in the presence of weak magnetic fields, differentially rotating gas disks subject to the 
robust magneto-rotational instability (MRI). This instability is caused by the back- reaction 
of magnetic fields to the fluid motion. Thus, in order to study the evolution of magnetic fields 
in differentially rotating disks, we should solve the momentum equation including Lorentz 
force simultaneously with the induction equation. This approach is called dynamical dy- 
namo. The developments of computational magnetohydrodynamics (MHD) enabled us to 
study this dynamical dynamo process by direct three-dimensional (3D) global simulations. 
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In the context of accretion disks, the nonlinear growth of MRI was first studied by local 
3D MHD simulations (e.g., Hawley, Gammie & Balbus 1995; Matsumoto & Tajima 1995; 
Brandenburg et al. 1995). They showed that MRI drives magnetic turbulence in accretion 
disks and amplifies magnetic fields. Subsequently, global 3D MHD simulations of initially 
weakly magnetized rotating torus were carried out (e.g., Matsumoto 1999; Hawley 2000; 
Machida, Hayashi & Matsumoto 2000). They showed that the magnetic field amplification 
saturates when f3 = P gas /Pmag ~ 10. The magnetic fields maintained in the quasi-steady 
state are mainly toroidal but have turbulent poloidal components. They also showed that 
in the innermost region of black hole accretion disks where radial infall of disk gas becomes 
important, spiral magnetic fields which accompany radial field reversals are created (e.g., 
Machida & Matsumoto 2003). 

In numerical studies of accretion disks, the gravitational field is assumed to be that of 
the central point mass. It is straightforward to extend this gravitational potential to more 
general potentials such as those of spiral galaxies. Since the galactic gaseous disk is also a 
differentially rotating disk, it subjects to the MRI. The growth of MRI and generation of 
MHD turbulence in galactic HI disks were discussed by Sellwood & Balbus (1999). Kim, 
Ostriker & Stone (2003) carried out three-dimensional local MHD simulations of galactic 
gaseous disks and showed that MRI drives interstellar turbulence and helps forming the 
giant molecular clouds. The coupling of the thermal instability and MRI in the interstellar 
space has been studied by Piontek & Ostriker (2004). 

Dziourkevitch & Elstner (2003) reported the results of global three-dimensional MHD 
simulations of the galactic gaseous disks including the dynamo a term. More recently, 
Dziourkevitch, Elstner & Riidiger (2004) carried out global three-dimensional MHD simula- 
tions of galactic gaseous disks without assuming the dynamo a term starting from vertical 
magnetic fields and showed that mean magnetic fields and interstellar turbulence are gener- 
ated. However, their simulation area was limited to 5kpc in the radial direction and lkpc in 
the vertical direction. Kitchatinov, & Riidiger (2004) carried out linear but global analysis 
for MRI in a disk geometry and showed that MRI can amplify very weak seed magnetic 
fields in proto-galaxies. 

In this paper, we report the results of global 3D MHD simulations of galactic gaseous 
disks by using the axisymmetric gravitational potential created by stars and the dark matter. 
We solve the nonlinear MHD equations without introducing the phenomenological dynamo 
a parameter. The initial magnetic field is assumed to be toroidal. This initial condition 
enables us to start the simulation without much disturbing the interface between the disk and 
the disk corona. When we start from vertical magnetic fields, efficient angular momentum 
loss near the surface of the initial torus drives an avalanche-like surface accretion, which 
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deforms magnetic field lines into an hourglass shape and produces magnetocentrifugally 
driven outflows (e.g., Matsumoto et al. 1996). For simplicity we neglect the multi-phase 
nature of galactic gaseous disk, non-axisymmetric spiral gravitational potential, self gravity 
and supernova explosions. These processes will be incorporated in subsequent papers. 

In section 2, we present numerical methods. Numerical results are given in section 3. 
Section 4 is devoted for summary and discussion. 



2. NUMERICAL METHODS 



2.1. Basic Equations 

The basic equations we solve are resistive MHD equations, 
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In these equations, p, p, B, v, <fi and 7 are the density, pressure, magnetic fields, velocity, 
gravitational potential and the specific heat ratio, respectively. The electric field E is related 
to the magnetic field B by E = (—vxB + 77V xB) jc. We assume the anomalous resistivity 

V 

r]o(v d /v c - l) 2 , (when v d > v c ), 
0, (whent> d < v c ), 

where v d = j/p is the electron-ion drift velocity and v c is the critical velocity above which 
anomalous resistivity sets in ( Yokoyama & Shibata 1994). 
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We adopted cylindrical coordinates (zu, ip, z). For normalization, we take the unit radius 
zu = 1 kpc, and unit velocity v = 255 km/s. The unit time is t — ^o/ v o = 3.8 x 10 6 years. 
For specific heat ratio and resistivity, we adopt 7 = 5/3 and i] = 0.01v G7 , respectively. 



2.2. Initial Condition 

As the initial density distribution, we assume a rotating equilibrium torus threaded by 
azimuthal magnetic fields embedded in the isothermal, non-rotating, spherical hot halo ( 
Okada, Fukue & Matsumoto 1989). We assume that the gas torus has constant specific 
angular momentum L and assume polytropic equation of state, P = Kp 1 where K is a 
constant. This initial condition simulates the situation that after a spiral galaxy is formed, 
the interstellar gas is supplied either by accretion of intergalactic gas which has some specific 
angular momentum or by supernova explosions following bursts of star formation inside the 
galaxy at some radius. As we shall show later, the angular momentum distribution of the 
gas settles into that determined by the galactic gravitational potential. Thus, the final 
distribution of density and angular momentum is not sensitive to the initial distribution. 

We assume that the Alfven speed va is a function of wB v , 

where B v is the toroidal magnetic field and 7i and \i are constants. We take \x = 7. Using 
these assumptions, we can integrate the equation of motion into the potential form, 

= constant, (7) 

where v s = ('jP/p) 1 ^ 2 is the sound speed and ^5 = ty(mb,0). We take the reference radius 
w\> at the radius of the density maximum of the torus. We adopt Wb = 10 kpc and set 
p{wbi 0) = Pb- Using equation (7), we obtain the density distribution 



P = Pb 



1 1/(7-1) 

max{* b - <f>(w, z) - L 2 /{2w 2 ), 0} 



.^{7/(7-1)} (i + /3 & " 1 ^- 1 )M 2(7 - 1) )_ 

where (3 b — {^K/l-Cj/wf 1 ' 1 ^ is the ratio of gas pressure to magnetic pressure at {m,z) 
(zub, 0). The thermal energy of the torus is parameterized by 



(8) 



E th = ^ = 0.05, (9) 

7^o 
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where c s & is the sound speed at (zu,z) = (zub,0). This sound speed corresponds to the 
temperature T5 ~ 2 x 10 5 K. This means that we assume a warm torus as the initial 
condition. For unit density, we take pb — 3 x 1CT 25 g/cm 3 . The initial gas pressure at 
(zu, z) = (zub, 0) is 

P h = P(zu b ,0) = - Pb c 2 sb 

7 

= 10- 11 ( Pb - - ) ( \-] dyne/cm 2 . (10) 

\3 x 10" 25 g/cm 3 y V 2 x 10 K / 



The density of the corona p c is given by 



10 4 pb exp 



c 



(11) 



where c sc is the sound speed in the corona and (pb is the gravitational potential at (zu, z) = 
(zub,0). We take the characteristic coronal pressure P c = p c c 2 sc ~ 10 _4 pf> c L = 3 x 10~ 3 P&. 

As the gravitational potential <fi, we adopted a model of our galaxy given by Miyamoto 
& Nagai (1975) and Miyamoto, Satoh & Ohashi (1980). They assumed centrifugal balance 
only in the galactic plane. To reproduce the galactic rotation curve, three different sets of 
parameters (a\, b±, Mi), (a 2 , b 2 , M 2 ), (a 3 , 63, M 3 ) are introduced in the gravitational potential 
as follows, 

3 GM 
^•^) = E [g ,2 + {oi + ( ,2 , + 6 y 5}T s - (12) 

Here i — 1 and 2 take into account the gravitational field by bulge and disk stars, respectively. 
The third component i = 3 is the contribution by dark matter. In this equation, the constants 
dj, bi and Mj have dimension of length and mass, respectively, and G is the gravitational 
constant. This potential is axisymmetric. The values of the constants are given in Table 1. 
Model I and II are models without dark matter and model III- VII include the dark matter. 

Figure 1 shows the isocontours of the gravitational potential for model I. Figure 2 shows 
the isocontours of the initial density and toroidal magnetic field for model I. 



2.3. Simulation Code 

We solved the MHD equations by using a three-dimensional MHD code based on a mod- 
ified Lax-Wendroff method ( Rubin & Burstein 1967) with artificial viscosity ( Richtmyer & 
Morton 1967). This code was originally developed to carry out two-dimensional axisymmet- 
ric global MHD simulations of jet formation from accretion disks (e.g., Shibata & Uchida 
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1985; Matsumoto et al. 1996; Hayashi, Shibata & Matsumoto 1996). The MHD code was 
extended to three- dimensions and has been used to carry out global three-dimensional MHD 
simulations of accretion disks (Machida et al. 2000; Machida & Matsumoto 2003; Kato, 
Mineshige & Shibata 2004). 



The number of mesh points is (N w , N v , N z ) = (250, 64, 319) for model I- IV. The grid size 
is Aw = 0.05, Az = 0.01 for < w < 6.0, < z < 2.0, and otherwise increases with w and z 
as Azk + i = 1.05Azfc. We used this nonuniform mesh to concentrate meshes near the galactic 
plane and near the central region. The size of the simulation box is kpc < w < 56 kpc, 
and kpc < z < 10 kpc. We used this large simulation box in order to simulate the buoyant 
rise of magnetic flux into the galactic corona and to reduce the effects of the outer boundary 
condition. 

We imposed free boundary conditions at w = 56 kpc and z — 10 kpc where waves can 
be transmitted. For azimuthal direction, we include full circle (0 < tp < 2n) for models 
I-IV. We carried out full circle simulations because we are interested in the formation of 
large-scale mean magnetic fields. We should note that the amplitude of modes with low 
azimuthal mode numbers (m = 1 or m = 2) are often large in global simulations of accretion 
disks (e.g., Machida & Matsumoto 2003). In order to study the effects of the azimuthal 
resolution and the azimuthal size of the simulation region, we carried out simulations for 
model V, VI and VII, in which the azimuthal simulation region is limited to < tp < tt/2. 
Periodic boundary conditions are imposed in the azimuthal direction. The other parameters 
are the same as those for model III (Table 2). 

We assume symmetric boundary condition at the equatorial plane z = and simulated 
only the upper half space. The physical quantities at rotation axis (w = 0) are computed 
by averaging those at the mesh point next to the axis. In models I and III- VII, we imposed 
absorbing boundary condition at r = [w 2 + z 2 ) 1 ^ 2 = r in = 0.8 kpc. We introduce a damping 
layer inside r = r in . In this layer, the deviation of physical quantity q from initial value qo 
is reduced after each time step with a damping rate a, 



This damping layer serves as the non-reflecting boundary which absorbs accreting mass and 
waves propagating into r = 0.8. This boundary condition is adopted in our standard model 



2.4. Simulation Region and Boundary Conditions 



q' = q - a (q- q ). 



(13) 



Here we take 




(14) 
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because the gas accreted to the central region of the galaxy will be converted to stars or 
absorbed by the central massive black hole. In model II, we imposed no specific condition 
at r = r in . Thus the accreting mass piles up in the central region. 

2.5. Initial Perturbation 

To initiate the evolution, we imposed random perturbation for the azimuthal velocity 
with maximum amplitude O.Olfo. We checked the stability of the initial state by carrying 
out numerical simulations without including magnetic fields and without imposing velocity 
perturbations. Since the boundary between the rotating disk and the static corona is not 
exactly in hydrostatic equilibrium, small motions with maximum velocity 0.05t>o appears in 
this region. The radial velocity in the dense region of the torus is less than O.OOlfo. The 
radial flow in the boundary layer between the disk and the corona slightly modifies the 
density distribution and smears the difference of the rotation speed between the disk and 
the corona. However, the flow does not significantly modify the density distribution of the 
main part of the torus. At t — 2 Gyr, the maximum radial speed is about 0.0045t>o even at 
the interface between the disk and the corona. The torus stays in the equilibrium state for 
time scale longer than 2.3 Gyrs. 

3. NUMERICAL RESULTS 

3.1. Time Evolution of the Model with Absorbing Inner Boundary Condition 

First we present the results of Model I (the standard model). The inner boundary is 
treated as an absorbing boundary where accreting mass and magnetic fields are absorbed. 
Figure 3 shows the time evolution of model I. Color contours show density distribution and 
solid curves show magnetic field lines. The left panels show the isosurface of the density 
(p/pb = 10~ L5 ), density distribution at the equatorial plane and magnetic field lines. We 
plotted the region —15 kpc < x,y,z < 15 kpc. The right panels show the equatorial density 
and magnetic field lines projected onto the equatorial plane. As the MRI grows, matter 
accretes to the central region by losing the angular momentum. The initial torus deforms 
itself into a flat disk. On the other hand, the matter in the outer part of the disk gets 
angular momentum and expands. The magnetic fields globally show spiral structure and 
locally show turbulent structure. 

Figure 4 shows the close up view of the density distribution and magnetic field lines 
projected onto the equatorial plane at t = 1000t (= 3.8 Gyr). The box size is —15 kpc < 
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x,y < 15 kpc. The density distribution shows weak non-axisymmetric spiral pattern. 

Figure 5a shows the time evolution of magnetic energy for model I averaged in 2 kpc < 
zu < 5 kpc, kpc < z < 1 kpc and < p < 2n. Figure 56 shows the plasma (3 averaged in 
the same region. As the MRI grows in the disk, the magnetic energy increases and its strength 
saturates when log[(S 2 /87r)/Pfe] ~ —1.5. For our Galaxy, this saturation level corresponds to 
B ~ 3/x G. The magnetic energy is maintained for time scale at least 2 Gyrs. As magnetic 
energy is amplified, the average value of plasma j3 = (P) / (B 2 /87r) decreases and stays around 
(f3) ~ 20 where (f3) is the spatial average of local plasma f3 in 2.0 kpc < w < 5.0 kpc, 
kpc < z < 1.0 kpc and < ip < 2n. This saturation level of magnetic energy is smaller 
than that observed in our galaxy. Further amplification of magnetic energy by, for example, 
supernova explosions may be necessary to explain the strength of magnetic fields in our 
galaxy. 

Magnetic field lines shown in the left panels of figure 3 indicate that magnetic fields 
initially confined in the disk emerge from the disk and form extended magnetized corona. 
Figure 6 shows the vertical distribution of the gas pressure and the magnetic pressure at 
zu = 3 kpc, 5 kpc and 10 kpc at t — and t = 1000t (= 3.8 Gyr). The magnetic flux 
initially confined in the torus fills almost the whole simulation region. The plasma j3 in the 
corona is 1 < f3 < 10 in 1 kpc < z < 10 kpc. 

Magnetic fields emerge from the disk to the corona due to the growth of the Parker 
instability (Parker 1966). Nonlinear growth of Parker instability in disks was studied by 
Matsumoto et al. (1988) by two-dimensional MHD simulations taking into account the ver- 
tical variation of the gravitational field. The effects of the differential rotation and MRI on 
the Parker instability were studied by Foglizzo & Tagger (1994, 1995). Miller & Stone 
(2000) carried out local three-dimensional MHD simulations of gravitationally stratified dif- 
ferentially rotating disks and showed that strongly magnetized corona is created due to the 
buoyant rise of magnetic flux from the disk to the corona. Machida et al. (2000) showed by 
global three-dimensional MHD simulations that buoyantly rising magnetic loops are formed 
above accretion disks. Our numerical results are consistent with these previous results. 

Figure 7 shows the radial distribution of azimuthal velocity v v and the density p averaged 
in kpc < z < 0.3 kpc and < p < 2n at t = 3.8 Gyr. In the inner region (zu < 10 kpc), 
the radial profile of azimuthal velocity of model I is similar to that observed in our Galaxy. 
The density distribution in Figure 7b indicates that the initial torus spreads and becomes 
flat in the inner region. 

Figure 8a shows the time evolution of the ratio of Maxwell stress averaged in 2 kpc < 
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tu < 5 kpc, kpc < z < 1 kpc and < ip < 2tt to the initial gas pressure at zu = 10 kpc, 

H-w>- (i5) 

The Maxwell stress increases as the gas accretes and later decreases gradually. 

Figure 8b shows the time evolution of the accretion rate at zu = 2.5kpc defined by 

p2ir p2.5kpc 

M^ =2 .5 = - / / pv m zudtpdz. (16) 
Jo Jo 

The unit of the accretion rate is 1M Q / yr when pi, = 3 x 10 -25 g/cm 3 . Numerical result 
indicates that M ~ 1O~ 3 M / yr. Since the mass of the initial torus is M torus = 5 x 10 8 M Q , 
mass accretion takes place for time scale longer than 10 Gyr. 



3.2. Effects of the Central Absorber and Dark Matter 

In this subsection, we show the effects of the central absorber and the dark matter. In 
model I, we assumed that the interstellar matter accreting to the central region are absorbed 
at r = 0.8 kpc. This mimics the existence of the central black hole or the efficient conversion 
of gas to stars in the central region. When such absorber does not exist, accreting matter 
piles up in the central region. As a reference model, we carried out a simulation of model II, 
in which we allow the accreting mass to pile up in the central region. 

Figure 9 shows the density distribution and magnetic field lines at t — 1000to(= 3.8 Gyr) 
for model II. Since accreted matter accumulates in the central region, it forms a high density 
bulge. 

Figure 10 shows the distribution of density and magnetic field lines at t — 3.8 Gyr for 
model III. In this model, the effect of dark matter is included. Absorbing boundary condition 
is imposed at r = 0.8 kpc. The density distribution and magnetic field lines are similar to 
those of model I. 

Figures 11a and 116 show the time evolution of magnetic energy and plasma (3 averaged 
in 2 kpc < zu < 5 kpc, kpc < z < 1 kpc and < tp < 2ir, for models I, II and III. The 
magnetic energy saturates at approximately the same level for all models. The magnetic 
energy is maintained more than 2Gyrs. However, after significant amount of gas accretes, 
the value of plasma (3 in model II becomes larger than that of other models, because in model 
II, the gas density and pressure increase with time. Figure 11c shows the ratio of Maxwell 
stress to gas pressure averaged in 2 kpc < zu < 5 kpc, kpc < z < 1 kpc and < <p < 2n 
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for model I, II and III. The time evolution of Maxwell stress has almost the same profile in 
all models. Figure lid shows the time evolution of the accretion rate at w — 2.5 kpc for 
these models. In model II, the accretion rate is about 10 times larger than that of model I 
and model III because the density near w = 2.5 kpc of model II is larger than the density 
of model I and model III (see figure 126). 

Figure 12 shows the distribution of azimuthal velocity and density p averaged in 

kpc < z < 0.3 kpc and < tp < 2tt at t = 3.8 Gyr. The azimuthal velocity of model 

1 and model III is similar to the observed profile of our Galaxy in the inner region. The 
velocity profile of model III coincides with the observed profile in all region because the 
effects of dark matter are included. On the other hand, the velocity distribution of model II 
is different from the observed rotation curve. In the inner region, the density of model II is 
ten times larger than that of model I and model III because the accreted gas piles up in the 
central region. When mass piles up in the central region, since the pressure gradient force 
in the radial direction reduces the effective gravity in the radial direction, the equilibrium 
rotation speed becomes much less than the observed rotation speed. Thus, the existence 
of the central absorber or conversion of infalling gas to stars are essential to reproduce the 
observed rotation curve. In the outer region (8 kpc < w) velocity and density of model I 
and model II almost coincide. 



3.3. Spatial and Temporal Reversals of Mean Magnetic Fields 

Figure 13 shows the magnetic field lines depicted by mean magnetic fields for model 
III at t — 3.8 Gyr. Regions colored in orange or blue show domains where mean azimuthal 
magnetic field at z — 0.25 kpc is positive or negative, respectively. The box size is 30 kpc x 
30 kpc. We computed the mean magnetic field at each grid point 

Ei+10 sr^3+ 3 v^ fe + 10 E>/ \ 
Z=i-10 l^m=j-3 Z^n=fc-10 - tf \ W h Vm, Z n ) 

<Pj, *k) = 21x7x21 > ( 17 ) 

by averaging magnetic fields inside ±10 meshes in w and z directions, and ±3 meshes in 
ip direction. The mean magnetic field is mainly toroidal but shows reversals of toroidal 
components around w ~ 5 kpc at this time. 

The azimuthal direction of mean magnetic fields in the disk changes with radius in 
regions where channel-like flow develops. Hawley & Balbus (1992) showed by axisymmetric 
simulations that channel-like flows appear in the nonlinear stage of MRI. Even when the 
unperturbed magnetic field is purely azimuthal, spiral channel flows appear as a result of 
the nonlinear growth of the non-axisymmetric MRI ( Machida & Matsumoto 2003). Since 
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such magnetic channels move inward as the infalling gas loses angular momentum, the radius 
of the field reversal changes with time. 

Figure 14 shows the azimuthally averaged mean toroidal field (B v ) (mean field) and the 
standard deviation of azimuthal field a/ ((B 9 — B^) 2 ) (fluctuating field) at t — 3.1 Gyr where 
( ) denotes the spatial average. The dashed, dotted and solid curves show the distribution 
of mean azimuthal magnetic fields averaged in azimuthal direction (0 < tp < 2tt) and in 
0.1 kpc < z < 1.0 kpc (equatorial region), 1.0 kpc < z < 3.0 kpc (coronal region) and 
0.1 kpc < z < 3.0 kpc, respectively. The vertical bars show the standard deviation of 
the azimuthal field. The strength of fluctuating field is comparable to the mean magnetic 
field. The dash-dotted curve shows the initial distribution of azimuthal mean magnetic fields 
averaged in 0.3kpc < z < 2.0kpc and < (p < 2tc. The strength of equatorial azimuthal 
magnetic field (0.1 kpc < z < 1 kpc) is larger than the initial azimuthal field by factor 2. 

The mean azimuthal magnetic field in the coronal region (1.0 kpc < z < 3.0 kpc) has 
direction opposite to that in the disk region. The total magnetic flux threading the zu-z 
plane is almost conserved. 

Figure 15 shows the spatial distribution of mean azimuthal magnetic fields for model III 
at t — 590to(= 2.2 Gyr) and 826to(= 3.1 Gyr). Blue and orange indicate regions where mean 
magnetic field B y threading the y = plane is positive or negative, respectively. Arrows 
show the direction of magnetic fields. These panels show that the mean azimuthal field B y 
threading the y = plane reverses its direction in the coronal region. The maximum height 
of the magnetized region (the wavefront of the azimuthal magnetic field) at w — 10 kpc 
locates at z = 4 kpc at t — 2.2 Gyr and z ~ 10 kpc at t — 3.1 Gyr. It indicates that the 
azimuthal magnetic flux is rising. 

The seeds for the reversals of the azimuthal magnetic fields are generated during the 
growth of MRI because when the initial magnetic field is purely azimuthal, the growth rate 
of the MRI is larger when q = (k y /k z ) 2 is small, where k y and k z are the wavenumber in 
the azimuthal direction and the vertical direction, respectively (e.g., Matsumoto & Tajima 
1995). Thus, the azimuthal direction of the perturbed magnetic field changes with height. 
In the nonlinear stage when the vertical magnetic fields are produced by the buoyant rise 
of the azimuthal magnetic flux, the formation of the channel flows (e.g., Hawley & Balbus 
1992) due to the nonlinear growth of MRI also produces the reversals of azimuthal magnetic 
fields. 

Figure 16a shows the time variation of the vertical distribution of azimuthally averaged 
magnetic field {B v ) for model III at zu — 10 kpc. The regions where {B v ) < are not 
plotted. The wavefront of the azimuthal flux propagates in the vertical direction with speed 
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v ~ 8 kpc/1.5 Gyr ~ 5 km/s. This speed is comparable to the Alfven speed. The vertical 
distribution of the azimuthal field (B^) approximately follows the self-similar solution of the 
nonlinear Parker instability (B v ) oc z' 1 (Shibata et al. 1989). Figure 166 shows the time 
evolution of the height of the wavefront of the azimuthal magnetic flux (solid curve) and 
the height of reversals of azimuthal magnetic fields (dashed curve and dash-dotted curve) 
at w = 10 kpc. The exponential increase of the height of the wavefront is consistent with 
the self-similar solutions of the nonlinear Parker instability (Shibata et al. 1989). After 
an azimuthal magnetic flux with one direction rises, another azimuthal magnetic flux with 
opposite direction rises. The MRI unstable disk continuously feeds the corona with magnetic 
fields. 

In figure 17 we plotted the time variation of mean azimuthal field of model III. Magnetic 
fields are averaged in azimuthal direction (0 < ip < 2n), radial direction (5 kpc < zu < 
6 kpc), and in 0.1 kpc < z < 1.0 kpc (equatorial region), 1.0 kpc < z < 3.0 kpc (coronal 
region) and 0.1 kpc < z < 3.0 kpc (disk + corona), respectively. The equatorial magnetic 
field in 5 kpc < w < 6 kpc changes direction with time. The azimuthal magnetic fields 
reverse their direction with timescale t ~ 300to(~ 1 Gyr). This timescale is comparable to 
that of the buoyant rise of azimuthal magnetic flux, which takes place in t ~ 10H/va, (e.g., 
Parker 1966), where H is the half thickness of the disk (~ 0.5 kpc), and va is the Alfven 
speed around the disk-corona interface (va ~ 5 km/s). 

Here we would like to point out that the amplification of equatorial magnetic field is 
enabled by the escape of magnetic flux from the disk to the corona. When the total azimuthal 
magnetic flux is conserved, the azimuthal magnetic flux inside the disk can increase when 
magnetic flux escapes from the disk to the corona. Although azimuthal magnetic flux is not 
exactly conserved in resistive MHD simulations, the flux is nearly conserved because r\ is 
small. Thus the buoyant escape of azimuthal magnetic flux from the disk to the corona is 
compensated by the amplification of azimuthal magnetic flux inside the disk with polarity 
opposite to that in the corona. 

3.4. Dependence on the Strength of Initial Magnetic Fields 

We also carried out a simulation (model IV) starting from the initially weaker magnetic 
field (f3 b = 1000). In this model, the magnetic energy saturates at almost the same level 
as in model III [B^ ~ l/i G at zu ~ 5 kpc). Figure 18 shows the spatial distribution of 
mean azimuthal magnetic fields for model IV at t = 824to(= 3.1 Gyr) and at t = 1115to(= 
4.2 Gyr). Blue and orange indicate regions where mean magnetic field B^ threading the 
y = plane is positive or negative, respectively. Again, the azimuthal magnetic field changes 
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direction with radius and with height. The striped distribution of azimuthal magnetic fields 
and the vertical propagation of the reversals of azimuthal magnetic fields are similar to those 
in model III (figure 15). 

In Figure 19 we plotted the time variation of mean azimuthal field of model IV. Magnetic 
fields are averaged in azimuthal direction (0 < ip < 2ir), radial direction (5 kpc < w < 
6 kpc), and in 0.1 kpc < z < 1.0 kpc (equatorial region), 1.0 kpc < z < 3.0 kpc (coronal 
region) and 0.1 kpc < z < 3.0 kpc (disk + corona), respectively. These results indicate 
that the direction of mean azimuthal magnetic field reverses several times. The amplitude 
of oscillation of azimuthal magnetic fields increases with time. The interval of reversal of 
azimuthal magnetic fields is ~ 1 Gyr (2 Gyr for a full cycle). This timescale coincides with 
that in model III, in which the initial plasma j3 is 100. The direction of azimuthal magnetic 
field in the coronal region is opposite to that of the equatorial region. 

3.5. Dependence on the Size of the Azimuthal Region and Mesh Size 

In order to study the dependence of numerical results on the size of the azimuthal 
simulation region and mesh sizes in the azimuthal direction, we carried out simulations in 
which we used 64 (model V), 32 (model VI) and 16 (model VII) grid points in < (p < n/2. 
Other parameters are the same as those in model III. 

Figure 20a shows the early stage (t < 600to) of the magnetic field amplification in 
2 kpc < w < 5 kpc and figure 20 b shows the time variation of magnetic energy in the non- 
linear saturation stage. We found that in the early stage (t < 600t ), the growth of the 
magnetic energy in model VII is similar to that of model III, which has the same azimuthal 
grid resolution. In the nonlinear saturation stage (figure 206), however, the magnetic energy 
in model VII decreases and becomes smaller than that in model III. These results are con- 
sistent with the results of local 3D MHD simulations of accretion disks ( Hawley, Gammie 
& Balbus 1995), in which they showed that the nonlinear saturation level of the magnetic 
energy increases with the size of the simulation region. 

When the azimuthal size of the simulation region is fixed (model V, VI and VII), the 
saturation level of the magnetic energy (t > 600to) is larger in the high resolution model 
(model V) than that in low resolution models (model VI and VII). This result is also consis- 
tent with the results of local 3D MHD simulations of the growth of MRI in accretion disks 
(e.g., Hawley, Gammie & Balbus 1995, 1996). 

The time evolution of the magnetic energy before its saturation (figure 20 a) depends 
on how the gas infalls and carries in magnetic fields. The earlier rise of magnetic energy in 
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the highest resolution model (model V) is due to the faster growth of the MRI in the torus. 
The growth rate of the MRI is sensitive to the azimuthal grid resolution because the most 
unstable wavenumber k» ~ Q/va, where k» is the wavenumber parallel to the unperturbed 
magnetic fields (e.g., Balbus & Hawley 1992) is hardly resolved in low resolution simulations. 
When (3 ~ 100, the most unstable wavelength is \ ma x ~ 2itva/Q ~ 2irH(2/ (3) 1 / 2 ~ H where 
H is the disk half thickness. Thus, \ max ~ 0.5 kpc ~ (tt/2) x 10 kpc/30. We need at least 
120 mesh points in0<y?<7r/2 when one wavelength is resolved by 4 mesh points. Thus, in 
all models we adopted in this paper, the fastest growing mode is not numerically resolved. 
These estimates indicate that the growth rate of MRI is affected by the grid size and that 
the growth rate in models III, VI and VII at zu — 10 kpc are smaller than that in model V. 
Therefore, the infall of the torus gas is delayed in low-resolution models. 

Let us discuss why magnetic energy in low resolution models (model III and VII) is 
larger than that in high resolution models (model V and VI) in the nonlinear stage before 
the saturation (t < 600to)- These magnetic fields are carried into 2 kpc < w < 5 kpc with 
the gas infalling along the channels of magnetic fields. Since non-axisymmetric parasitic 
instabilities ( Goodman & Xu 1994) and magnetic reconnection (e.g., Sano & Inutsuka 
2001; Machida & Matsumoto 2003) which break up the channel flow can be captured better 
in high resolution simulations, the growth rate of the magnetic energy becomes smaller in 
high resolution models than that in low resolution models. 

Figure 21a compares the vertical distribution of azimuthal magnetic fields at w — 10 kpc 
at t — 800to in model III, V, VI and VII. The region where B v < is not plotted. The 
azimuthal magnetic flux buoyantly rises into the corona and forms extended magnetized 
region. Numerical results indicate that the magnetic flux rises faster into the corona in 
the high resolution model (model V) than low resolution models (model VI and VII). High 
resolution simulations could resolve the most unstable wavelength of the Parker instability, 
A ~ IQH ~ 5 kpc in our simulation model (e.g., Parker 1966; Matsumoto et al. 1988). Thus, 
one wavelength of the Parker unstable loop is well resolved in high resolution model with 
64 grid points in0<(/9<7r/2attu = 10 kpc (model V) but hardly resolved in the low 
resolution model (model VII). This may be the reason why toroidal magnetic flux rises faster 
in model V than in other models (figure 216). 

When the azimuthal resolution is fixed (model III and VII), the magnetic flux rises 
faster in the model with larger azimuthal simulation region (model III). This dependence on 
the azimuthal simulation region comes from the horizontal expansion of the rising magnetic 
loops (e.g., Shibata et al. 1989). When the horizontal simulation region is restricted, since 
magnetic tension prevents the escape of the magnetic flux, the speed of the buoyantly rising 
magnetic loops is reduced. 
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3.6. Difference Between the Rotation Curve of Dark Matter and Gas 

Figure 22 compares the rotation curve of stars and dark matter computed from the 
gravitational potential (dashed curve) and the rotation curves for gas (solid curve) obtained 
from simulations for model III. The rotation curves of our Galaxy interior to the Sun has 
been obtained by HI gas (e.g., Gunn, Knapp & Tremaine 1979), and CO ( Burton & 
Gordon 1978) observations. The rotation speed in 5 — 10 kpc from the galactic center is 
220 — 260 km/s and similar to the rotation speed shown in figure 22. The rotation speed of 
our galaxy distant from the galactic center [w > 10 kpc) is obtained by radio observation of 
HI gas (e.g., Merrifield 1992). The rotation speed is nearly flat and v v ~ 240 km/s. The flat 
rotation curve is similar to the rotation curve in the outer part of external spiral galaxies 
(e.g., Rubin et al. 1985). Our numerical result reproduces these flat rotation curves when 
we use the gravitational potential which includes the dark matter. 

Since the rotation curve of our galaxy distant from the galactic center is obtained by 
radio observations of HI gas, the observationally measured rotation curve is that for the gas. 
The rotation curve of the neutral hydrogen does not necessarily coincide with that of stars 
and the dark matter. Since they are almost frozen to the magnetic fields, if the magnetic 
arm rotates with angular speed different from the dark matter, the neutral hydrogen is not a 
good tracer of mass distribution in our Galaxy. Battaner et al. (1992) suggested that when 
magnetic pinch force is large enough, they can explain the observed rotation curve of the 
Galaxy without dark matter. On the other hand, Sanchez-Salcedo & Reyes-Ruiz (2004) 
reported that the contribution of the magnetic field to the rotation curve is up to 20 km/s 
and it is not sufficient to banish the effect of dark matter. Our numerical results indicate that 
even in the presence of dynamically non-negligible magnetic fields, gas component shows a 
flat rotation curve similar to that observed by neutral hydrogen. Numerical results for model 
I indicates that if the dark matter does not exist, gas rotation curve decreases with radius. 
Thus, the flat rotation curve obtained by observations of neutral hydrogen really indicates 
the existence of dark matter. 



4. Summary and Discussion 

We carried out global three-dimensional magnetohydrodynamic (MHD) simulations of 
galactic gaseous disks to investigate how the galactic magnetic fields are amplified and main- 
tained. We showed the dependence of numerical results on physical conditions in the inner- 
most region, gravitational potential, and strength of initial magnetic fields. 

In all models, the magnetorotational instability (MRI) grows in the disk and the mag- 
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netic energy is amplified and maintained for more than 5 Gyrs. The saturation level of 
magnetic fields, however, is smaller than that expected from observations (B ~ several ixG) 
at w — 10 kpc. The multiphase nature of interstellar space, supernova explosions, cosmic 
rays and/or non-axisymmetric spiral potential may be important to further amplify magnetic 
fields. 

Numerical results indicate that mean azimuthal field inside the galactic disk is ampli- 
fied and that its direction changes quasi-periodically. The amplitude of the oscillation of the 
azimuthal field increases with time. The interval of reversal of equatorial azimuthal field is 
~ 1 Gyr. This timescale is comparable to the timescale of the buoyant rise of azimuthal 
magnetic flux from the disk to the corona. The amplitude of the fluctuating field is compa- 
rable to the mean magnetic field. The strength of the mean field only weakly depends on 
the initial strength of magnetic fields. In the coronal region, the mean azimuthal field has 
direction opposite to the mean field of the disk. 

Mean magnetic fields obtained from numerical simulations also show reversals in the 
radial direction near the equatorial plane around w = 5 kpc at t — 3.8 Gyr. The radius 
of the field reversal changes with time because the spiral magnetic channels which produce 
the reversal of azimuthal fields in the radial direction move inward as the disk gas accretes 
toward the galactic center. 

In the Galactic plane, magnetic field reversals are observed by Rotation Measure (RMs) 
of pulsars (e.g., Han et al. 2002). The reversals of mean magnetic fields are observed in our 
Galaxy between the local Orion arm and the inner Sagittarius arm ( Simard-Normandin & 
Kronberg 1979). More recent analysis imply axisymmetric field with two reversals ( Rand & 
Kulkarni 1989; Rand & Lyne 1994; see also Beck et al. 1996). Simulation results of model 
III (Figure 13) show reversals of mean magnetic fields near the equatorial plane around 
w — 5 kpc. Note that RMs measure the mean magnetic field in some specific direction from 
the Earth. When we draw lines from a point in Figure 13, we can see field reversals in several 
directions. Since the magnetic fields inside the galactic disk is turbulent, more and more 
field reversals will be observed as angular resolution of RMs increase. 

In conventional theories of aw-dynamos, various kinds of spatial and temporal oscilla- 
tions appear depending on the dynamo number D = C a Cn where Cn = H 2 Q/r] (H is the 
thickness of the disk, Q is the rotation angular speed and rj is the magnetic diffusivity) and 
C a = Ha/ 7] (e.g., Stepinski & Levy 1988). In our global MHD simulations, we do not need 
to introduce the dynamo a parameter. The amplification of magnetic fields due to MRI and 
the reproduction of toroidal fields through gas motions are self-consistently incorporated in 
our simulations. We introduced anomalous resistivity rj to handle magnetic reconnection 
which takes place in local region where current sheets are formed, but the turbulent diffu- 
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sivity is not explicitly assumed. We showed by direct 3D simulations that mean toroidal 
magnetic fields show reversals both in space and time. 

Tout & Pringle (1992) constructed phenomenological models of disk dynamos by taking 
into account the growth of magnetic fields due to MRI, escape of magnetic flux by the Parker 
instability, and the dissipation of magnetic fields by magnetic reconnection. They showed 
that the magnetic fields in the disk oscillate quasi-periodically with period comparable to 
the time scale of the buoyant escape of the magnetic flux. This timescale is consistent with 
the timescale of field reversals obtained from our MHD simulations. 

Let us discuss the similarities and dissimilarities of the magnetic activity of galactic 
gas disks produced by our numerical simulations and those in the Sun. Both objects are 
differentially rotating plasma with high conductivity in which magnetic fields are almost 
frozen to the rotating plasma. Thus the magnetic fields can be deformed and amplified by 
the plasma motion. 

The Sun is a slow rotator in which the centrifugal force is much smaller than the gravity, 
thus its atmosphere is confined by the gravity. The convective motions and/or the differential 
rotation in the convection region may be the source of magnetic field amplification in the 
Sun. However, since the magnetic flux tube in the convection region buoyantly rises with 
a timescale of the order of a month (e. g., Parker 1975), much shorter than the solar 
magnetic cycle (~ 22 year), solar dynamo activity is considered to be operating in the 
interface between the convection zone and the radiative zone, where magnetic buoyancy is 
small and differential rotation exists (e. g., Spiegel & Weiss 1980; Schmitt & Rosner 
1983). On the other hand, galactic gas disks are supported by the rotation and subject to 
differential rotation, which induces the growth of MRI in the whole disk. In galactic gas 
disks, since the gas disk is convectively stable, the buoyant escape of the magnetic flux takes 
place with timescale tb uoy ~ 10H/va ~ 10(l/S7)(/5 /2) 1 / 2 . This timescale is longer than the 
growth time of MRI turn ~ 1/^- Thus the magnetic fields can be amplified in the disk 
before escaping to the corona. 

Numerical results indicate that the total azimuthal magnetic flux is nearly conserved. 
When turbulent diffusivity exists, the azimuthal magnetic flux is not necessarily conserved. 
In our simulations, however, since magnetic diffusion only exists in the local current sheet, 
the magnetic flux is almost conserved. We found that the toroidal magnetic flux in the 
disk buoyantly rises into the corona. This result indicates that the buoyant escape of the 
magnetic flux from the disk enables the amplification of mean magnetic fields in the galactic 
gas disk. When the magnetic diffusion is not negligible, the total azimuthal magnetic flux 
is not conserved but the magnetic helicity is conserved. Blackman & Brandenburg (2003) 
pointed out the importance of helicity conservation and the role of coronal mass ejections, 
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which help sustaining the solar dynamo cycle. In our simulations, buoyantly rising magnetic 
loops carry helicity as well as the magnetic flux threading the disk. This escape of magnetic 
flux and helicity enable the disk to amplify the azimuthal magnetic flux in the equatorial 
region. 

We showed that when we remove the central absorber at r = 0.8 kpc, formation of 
dense gas bulge forces the rotation curve to deviate from that expected from the gravitational 
potential (model II). This indicates that in the central region of the galaxy, the gas absorption 
or the conversion of the accumulated gas to stars are essential to explain the observed rotation 
curve. The mass accretion rate to the central region is ~ 10~ 3 M Q / yr when we adopt the 
absorbing boundary condition. This accretion rate depends on the initial density of the torus, 
which can be formed by the infall of intergalactic matter or by the supernova explosions after 
bursts of star formation in the certain radius of the disk. When the galaxy was more gas 
rich in the early stage of its evolution, the accretion rate should be much higher. 

We also showed that the gas rotation curve approximately coincides with the rotation 
curve for stars and dark matter even when magnetic fields are dynamically important. This 
justifies us to use the gas rotation curve to estimate the distribution of the dark matter. 
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IAU XXXth General Assembly, 2003 and the workshop on Magnetic Fields in the Universe, 
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priority research project of the Graduate School of Science and Technology, Chiba University, 
ACT-JST of the Japan Science and Technology corporation, and UK-Japan collaboration 
on magnetic activities of the Sun, Stars and accretion disks (P.I., K. Shibata and N. Weiss). 
Numerical computations were carried out on VPP5000 at the Astronomical Data Analysis 
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Table I. Model parameters and inner boundary conditions 



Models 


ai 


ai 


a 3 


h 


b 2 


h 


Mi 


M 2 


M 3 




Boundary 


Model I 


0.0 


7.258 


0.0 


0.495 


0.520 


0.0 


2.05 


25.47 


0.0 


100 


Absorbing 


Model II 


0.0 


7.258 


0.0 


0.495 


0.520 


0.0 


2.05 


25.47 


0.0 


100 


Non- Absorbing 


Model III 


0.0 


6.2 


0.0 


0.47 


0.15 


31.2 


1.95 


17.4 


73.5 


100 


Absorbing 


Model IV 


0.0 


6.2 


0.0 


0.47 


0.15 


31.2 


1.95 


17.4 


73.5 


1000 


Absorbing 


Model V-VII 


0.0 


6.2 


0.0 


0.47 


0.15 


31.2 


1.95 


17.4 


73.5 


100 


Absorbing 



Table 2. Azimuthal grid size, number of grid points in azimuthal direction (iV^), and the 



azimuthal size of the simulation region 



Models 




N v 


simulation region 


Model I-IV 


tt/32 


64 


< up < 2tt 


Model V 


tt/128 


64 


< ip < tt/2 


Model VI 


tt/64 


32 


< up < tt/2 


Model VII 


tt/32 


16 


< ip < tt/2 



Note. - Other parameters of model V, VI 
and VII are the same as those in model III. 
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Galactocentric distance (kpc) 



Fig. 1. — Isocontours of the gravitational potential in Miyamoto & Nagai's (1975) model. 
This potential for model I does not include the effects of dark matter. 
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Fig. 2. — Isocontours of initial density distribution (upper panel) and azimuthal magnetic 
field (bottom panel) for model I. 
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Fig. 3. — Time evolution of density distribution log(p/p&) (color scale) and magnetic field 
lines (gray curves) in model I. The box size of left and right panels are 30 kpcx30 kpcx30 kpc 
and 30 kpc x 30 kpc, respectively. The right panels show the equatorial density and magnetic 
field lines projected onto the equatorial plane. 



Fig. 4. — Upper panel: The distribution of density log(p/p&) for model I in the equatorial 
plane at t = 3.8 Gyr. Lower panel: magnetic field lines in model I projected onto the 
equatorial plane at t — 3.8 Gyr. The box size is 30 kpc x 30 kpc. 
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Fig. 5. — (a)Time evolution of the magnetic energy and (b)volume averaged plasma (3(= 
(P) / (B 2 / 8n}) for model I averaged in 2 kpc < zu < 5 kpc, kpc < z < 1 kpc and < ip < 
2tt. 
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Fig. 6. — Vertical distribution of the gas pressure and magnetic pressure, (a) Initial distribu- 
tion of gas pressure, (b) Initial distribution of magnetic pressure (c) vertical distribution of 
gas pressure at t — 1000t and (d) the vertical distribution of magnetic pressure at t — 1000t - 
The solid, dashed and dotted curves indicate quantities averaged at w = 3 kpc, 5 kpc and 
10 kpc, respectively. 
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Fig. 7. — The distribution of (a)azimuthal velocity and (b)density p at t — 3.8 Gyr averaged 
in kpc < z < 0.3 kpc and < (p < 2n for model I. 




Fig. 8. — (a)Time evolution of Maxwell stress normalized by for model I. The Maxwell 
stress is averaged in 2 kpc < zu < 5 kpc, kpc < z < 1 kpc and < (p < 2n. (b)Time 
evolution of the accretion rate at 2.5 kpc from the center. The unit of the accretion rate is 
M = IMq/ yr when p b = 3 x 10~ 25 g/cm 3 . 
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Fig. 9. — The distribution of density (color scale) and magnetic field lines (gray curves) at 
t = 3.8 Gyr for model II. The box size is 30 kpc x 30 kpc. 
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Fig. 10. — The distribution of density (color scale) and magnetic field lines (gray curves) at 
t = 3.8 Gyr for model III. The box size is 30 kpc x 30 kpc. 
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Fig. 11. — Comparison of numerical results for model I, model II and model III. (a)Time 
evolution of the magnetic energy. (b)Time evolution of the plasma f3 averaged in 2 kpc < 
zu < 5 kpc, kpc < z < 1 kpc and < (p < 2n. (c)Maxwell stress averaged in the same 
region. The solid, dashed and dotted curves show the results for model I, model II and model 
III, respectively. (d)The time evolution of the accretion rate at 2.5 kpc from the center. 




Fig. 12. — (a)The radial distribution of azimuthal velocity v v and (b)radial distribution of 
density p averaged in kpc < z < 0.3 kpc and < if < 2n at t = 3.8 Gyr. The solid, dashed 
and dotted curves show the results for model I, model II and model III, respectively. 
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Fig. 13. — Magnetic field lines depicted by mean magnetic fields (solid curves) and isocon- 
tours of azimuthal component of magnetic fields at z = 0.25 kpc (color) at t = 3.8 Gyr for 
model III. Regions colored in orange or blue show domains where azimuthal magnetic field 
at z = 0.25 kpc is positive or negative, respectively. The box size is 30 kpc x 30 kpc. 
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Fig. 14. — Mean azimuthal field and its standard deviation (vertical bars) for model III at t — 
3.1 Gyr averaged in azimuthal direction (0 < ip < 2tt) and in 0.1 kpc < z < 1 kpc (dashed 
curve), 1.0 kpc < z < 3.0 kpc (dotted curve) and 0.1 kpc < z < 3.0 kpc (solid curve), 
respectively. The unit of is a/ Pbvl = 15// G in this model when p b = 3 x 10~ 25 g/cm 3 . 
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Fig. 15. — The distribution of mean azimuthal magnetic fields for model III at t — 590to(= 
2.2 Gyr) and 826io(= 3.1 Gyr). Blue and orange indicate regions where the mean field B y 
threading the y = plane is positive or negative, respectively. Arrows show the direction 
of mean magnetic fields. The mean azimuthal fields change their direction with height from 
the equatorial plane. 
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Fig. 16. — (a) Vertical distribution of azimuthally averaged magnetic field (B^) at w — 
10 kpc at t — 400to, 600to and 800to for model III. The dashed curve shows the initial profile 
of (Bp). The dotted curve shows the B v oc z^ 1 relation expected from the nonlinear theory of 
the Parker instability, (b) The time evolution of the wave front of rising magnetic flux where 
B v = 0.003a/ Pfct>o (solid curve) and the height of where direction of azimuthal magnetic 
field reverses (dashed curve and dash-dotted curve) for model III. The dotted curve shows 
the relation z — 1 + exp {(t — 300) /200]} expected from the nonlinear theory of the Parker 
instability. 
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Fig. 17. — The time variation of mean azimuthal field of model III averaged in 5 kpc < 
zxj < 6 kpc, < ip < 27T and 0.1 kpc < z < 1 kpc (dashed curve), 1.0 kpc < z < 3.0 kpc 
(dotted curve) and 0.1 kpc < z < 3.0 kpc (solid curve), respectively. The unit of is 
\J pbVQ = 15yU G in this model when pi, = 3 x 10~ 25 g/cm 3 . 
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Fig. 18. — The distribution of mean azimuthal magnetic fields for model IV = 1000) at 
(a) t = 824to(= 3.1 Gyr) and (b) t = 1115£o(= 4.2 Gyr). Blue and orange indicate regions 
where the mean field B y threading the y = plane is positive or negative, respectively. 
Arrows show the direction of mean magnetic fields. The mean azimuthal fields change their 
direction with height from the equatorial plane. 
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Fig. 19. — The time variation of mean azimuthal field of model IV averaged in 5 kpc < 
w < 6 kpc, < tp < 27T and 0.1 kpc < z < 1 kpc (dashed curve), 1.0 kpc < z < 3.0 kpc 
(dotted curve) and 0.1 kpc < z < 3.0 kpc (solid curve), respectively. The unit of is 
\J pbVQ = 15yU G in this model when pi, = 3 x 10~ 25 g/cm 3 . 
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Fig. 20. — The dependence of azimuthally averaged magnetic energy on the azimuthal grid 
resolution and the size of the azimuthal simulation region, (a) Field amplification stage and 
(b) non-linear saturation stage. Magnetic energy is averaged in 2 kpc < w < 5 kpc, and 
kpc < z < 1 kpc. Azimuthal simulation region is 2ir for model III and ir/2 for model V, 
VI and VII. The azimuthal grid size is 7r/32 for model III, and 7r/128, tt/64 and 7r/32 for 
model V, VI and VII, respectively. 




Fig. 21. — (a) Vertical distribution of azimuthal magnetic field B^ at w = 10 kpc when 
t = 800to for model III, V, VI and VII. The dashed curve shows the initial profile of B v . (b) 
The time evolution of the rising magnetic flux for model III, V, VI and VII. The height of 
the wave front of the rising magnetic flux where B v = 0.003 a/p&^o i s plotted. 
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Fig. 22. — The rotation curve of stars and dark matter computed from the gravitational 
potential (dashed curve) and the rotation curve for gas (solid curve) obtained by simulations 
for model III at t = 3.8 Gyrs. 



